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Abstract 

A rational interpolation method for approximating a frequency response is presented. 
The method is based on a product formulation of finite differences, thereby avoiding the 
numerical problems incurred by near-equal-valued subtraction. Also, resonant pole and 
zero cancellation schemes are developed that increase the accuracy and efficiency of the 
interpolation method. Selection techniques of interpolation points are also discussed. 


1 Introduction 

Consider the linear time-invariant system given by the state-space model 

x = Ax + Bu (1) 

y = Cx (2) 

where A € ]R n * xn *, B € Et njtXnB , C € R ncXnx , and the state vector, input vector, and 
output vector, x, u, and y, respectively, are properly dimensioned. We shall refer to the 
matrices, A, B, and C, as the state coupling matrix, the input coupling matrix, and the 
output coupling matrix, respectively. 

The frequency response of such a modeled system is defined as the Laplace transform of 
the input-output relationship evaluated along the jw-axis, 

G(u) = C(jul - A)~ 1 B (3) 


where 

0 < u> < oo. 

♦This research was supported in part by the Air Force Office of Scientific Research under Contract No. 
AFOSR91-0240. 
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In this paper, a fast and reliable interpolation method to compute frequency response is 
presented. The basic idea of this method is based on the simple Taylor series approximation 

G(u + h) = T 0 + T 1 h + >-- + T k h k + E k (4) 

but considered in the general interpolation form with fc+l interpolation points ho, hi,..., h k : 

G(w + h) = Go + G\(h — ho) + • • • + G k (h — ho)(h — h\) — (A — h k -i) + E k . (5) 

The coefficient matrices, Go, G \, . . . , G k are of size nc x ns as is the truncation error E k . 
Therefore, the cost of evaluating the matrix polynomial approximation 


P k (h) = Go + G\(h — ho) H h G k (h — ho)(h — h\) ■ ■ *(A — hjt_i) (6) 

is just knBnc floating-point operations (flops). The cost of computing each coefficient ma- 
trix is approximately the same as evaluating G by the method that would normally be 
preferred. 

The polynomial interpolation scheme works well as long as u is not near a resonant pole or 
zero of the system. In order to avoid this problem, we introduce methods of preliminary 
pole and zero cancellation. These greatly increase the accuracy of the interpolation scheme 
while causing only a negligible increase in the cost of computing the coefficient matrices. 

We shall also discuss the implementation of this algorithm including ideas on the selection 
of interpolation points. 

2 Existing Frequency Response Methods 

2.1 Straightforward Computation 

An obvious method for computing frequency response for a system modeled in state-space 
form is first to perform an LU decomposition in order to solve the linear system 

(ju>I-A)X = B, (7) 

followed by a matrix multiplication involving the solution to (7), 

G(w) = CX. 

This method does not exploit any special structure, e.g., sparse or banded, and therefore 
would only be used for general systems. To compute a frequency response implementing 
this method, for just one value of u, approximately + j(ns + nc)n\ + nonane flops 
axe required. As the number of desired frequency points becomes large, the calculation of 
the entire frequency response becomes computationally intensive. 
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2.2 The Principal Vector Algorithm 

In order to reduce computation cost, several methods have been developed in order to re- 
duce the cost of solving the linear system (7) either by exploiting the structure of the state 
coupling matrix or by implementing a similarity transformation to put the matrix A into 
an exploitable form. One method of the latter variety is the Principal Vector Algorithm 
(PVA) [10]. 

The idea of the PVA is to initially transform the state coupling matrix into a Jordan 
Canonical Form (JCF). The algorithm uses the principal vectors to compute the JCF in a 
more accurate way than previous such algorithms. Let 

A = M~ l JM (8) 

where J is in Jordan form. If we substitute this identity into (3), the frequency response 
becomes 


G(u) = CMM~ l (juI - A)~ l MM~ l B 

= C(jul — J)~ X B. (9) 

The initial transformation using the PVA to compute the JCF requires only 0(n flops if 
the state coupling matrix is not defective while 0(n\) flops are required if the matrix A 
is defective. Note that this transformation only occurs once, thus the cost is only incurred 
once. The advantage occurs in computations at each frequency point where the cost is 
reduced to 0(nx + n^nsnc) flops in the nondefective case and 0{\n^ + TMUfinc) A°P S 
in the defective case. So the computational saving occurs after the computation of one 
frequency point in the former case and n>i frequency points in the latter case. 

Although this algorithm produces significant savings in the computational cost of a fre- 
quency response, it can also frequently encounter numerical instabilities. First, the JCF is 
extremely unstable. The slightest perturbation can change a defective matrix into a non- 
defective matrix. Another problem is that the similarity transform may be ill-conditioned 
with respect to inversion depending on the basis of eigenvectors. If they are guaranteed 
to form a matrix which is well- conditioned with respect to inversion as would occur if the 
matrix were normal, the algorithm is very effective. 

2.3 The Hessenberg Method 

Another algorithm which uses similarity transformations to put the state coupling ma- 
trix into an exploitable form is the Hessenberg Method [4]. This algorithm is the current 
standard for computing the frequency response for generic dense systems. The Hessenberg 
Method, as its name implies, performs an initial transformation on the state coupling matrix 
to reduce it to upper Hessenberg form. So in this case, we use the identity 

A = Q~ 1 HQ f 
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where H is in upper Hessenberg form, instead of the JCF identity (8), in the frequency 
response (3). ; 

As with PVA, this initial transformation is performed only once at the start of the algo- 
rithm at a cost of O(n^) flops. When this transformation is used, the cost of computing 
the frequency response at each value of w becomes 0(n^(ng -f 1) + ngnc) flops. Usually, 
ng < riyt so a significant reduction in computation can be realized. 

Fortunately, there always exists an orthogonal transformation to reduce the state coupling 
matrix into an upper Hessenberg form. This prevents ill-conditioning from being introduced 
into the calculations by the similarity transformation as can occur with the Principal Vector 
Algorithm. 

2.4 Sparse Systems 

Many of today’s large ordered systems are sparse systems. A sparse system is one whose 
modeling matrices have relatively few nonzero entries when compared to the total number 
of entries. In such cases the Hessenberg Method should not be used. Instead of maintaining 
sparsity, the initial transformation will create a large dense system which then must be 
solved. There exist many storage techniques for sparse matrices which require a significantly 
smaller amount of memory allocation than a full matrix of the same order would require. 
Also, sparse matrix algorithms have been developed to exploit sparsity in order to reduce 
the computational costs in comparison to their dense counterparts. (See [6], [9], and [7].) 
These algorithms attempt to prevent the cost of solving the linear system (7) from growing 
to O(n^) flops. 

2.5 Frequency Selection Routines 

The cost of computing an entire frequency response can also be reduced by eliminating 
needless recalculations or overcalculations in attempts to get a desired resolution in the 
solution. When the frequency mesh is too coarse to give the required information, usually 
the user recomputes the entire frequency response. Often, the response from the previously 
computed frequency values either is recalculated or just ignored in the new calculation. 
Also, many times the user creates a fine frequency point mesh across the entire frequency 
range. Usually, only in small subregions is the finer mesh needed. A coarser mesh would 
suffice over the rest of the frequency range. 

In an effort to eliminate these unnecessary calculations but still give the required accuracy, 
so-called adaptive routines have been developed. These routines adapt the frequency point s 
selection to the characteristics of the system being analyzed. 

One such adaptive scheme is similar in nature to the QUANC8 adaptive integration routine 
[1]. The basic idea is first to select the endpoints of an interval in the desired frequency 
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region. Then the frequency responses of the two points are compared. If the difference 
between their magnitudes or their phases is greater than specified tolerances, the interval 
is divided in half. Then the three points are compared. If their differences are outside 
the tolerances, the subintervals are again halved. This subinterval halving continues until 
the tolerances are met across the entire interval or until a specified number of frequency 
points has been calculated. A single-input single-output variation of a method based on 
subinterval halving has been implemented commercially [5]. 

The use of a priori information, e.g., the locations of poles and zeros of a system, can also 
be used in the choice of frequency locations. More points are placed in the areas where the 
poles and zeros of a given system have an effect. Fewer points are placed outside these areas. 
Such a method is now being implemented in a linear system package [2] to automatically 
choose the frequency range over which the frequency response is computed as well as to 
determine the number of points needed to be calculated. 

These adaptive schemes also can be combined to form hybrid routines. This would permit 
an initial placement of points with the a priori method and then create the frequency mesh 
to join the regions between the areas of the initial placement. 


3 Polynomial Interpolation 

In order to compute the coefficient matrices, G \ , . . . , G k , of the interpolation equation 

P*(h) = Go + Gi(h - h 0 ) + • • • + G k (h - h 0 )(h - hi) • • -(h - h k . i) (10) 


finite differences will be employed. The first-order difference is defined as 

, /ri i i Af(hi) - M(h 0 ) 



while higher-order differences are defined as 


A/[/io, hi, ..., h n ] — 


Af [hi, ..., h^] M[ho, ..., h„ — i] 

h n - h 0 


If we let 
where 


M(h) = (jhl-A 0 ) 


-l 


A 0 = -jwl + A, 

the k th -order interpolation approximation can be written as 
P k (h ) = C(M(ho) + Af[h 0 ,hi](h-ho) + --- 

+M[h 0 , hj, .... h*](h - ho)(h - hi) • • -(h - h*_i))£ 


( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 
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with the interpolation error 


E k = G(u + h) - P k (h) 

k 

= C(M[h 0 , hi, h k , h) IJ(h - hi))B. (16) 

{=0 

Now, for convenience, define 

-Pfc(h) = M(ho) + M[ho,hi](h — ho) + • •• 

+M[h 0 , hi,..., h k ](h - ho)(h - hi) ■ • .(h - h k . i) (17) 

and 

E k = M(h)-P k (h) 

k 

= M[h 0 ,hi,...,h fc ,h]n(*-M- ( 18 ) 

i=0 

Although finite differences have a certain elegance to their formulation, they can encounter 
numerical inaccuracies due to the subtraction of near-equal- valued quantities. An extreme 
example of this is the case in which all of the interpolation points are the same. In theory, 
the first-order difference is exactly the first derivative of M, but numerically it is useless. 
Fortunately, the differences of the resolvent function (13), can be expressed in matrix prod- 
uct forms which avoid these cancellation problems as the following theorem shows. 


Theorem 1 For the resolvent function, the matrix difference functions in (12) satisfy 

M[ho, hi, . . ., h m ] = ( -j) m M{ho)M(hi ) ■ • -M(h m ). (19) 

Proof. Using (13) 

M(hi)-M(h 0 ) = (jhil - A 0 )~ l - (jhol - A 0 )- 1 

= ( jhol - Ao)" 1 {jh 0 I -Ao- (jhil - A 0 )) (jhj - Ao)' 1 
= (- j)(h ! - h 0 )M(h 0 )M(hi). 

Thus the first finite difference becomes 

M[ho,hi] = -jM(ho)M(hi) 

which proves (19) for m = 1. Now suppose that (19) is true for m-1. Since M(h 0 ), . . . ,M(h m ) 
all commute with each other, we find that 

M[h 0 , hi , . . . , hm] 
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= (Af [hi, h m ] - M[ho, h m ~i])/(h m - ho) 

= ( -j)”- 1 (M(M • • • M(h m ) - M(ho) ■ • ’ M(h m _j)) /(ft m - h 0 ) 

= (-jT ' 1 (M(M) • • -Af(fc m _i)) (M(Am) - M(M) /(fcm - M 

= (-j) m (M(^) • • • M(h m _ x )) M(fto)M(fc m ) 

= (-irM(MM(A x )---M(A m ). 

Thus (19) is true for m and thus, by induction, the theorem is true. □ 

If we now substitute the resolvent identity (19) into (17) and (18) and use the commutative 
property of the resolvent functions, the interpolation approximation becomes 

P k (h) = M(h 0 ) + (-j)M(h 1 )M(h 0 )(h-ho) + ••• 

+(-j)*M(h*)M(h*_ x ) • -M(ho){h - ho) ■••(/» - hk-i) (20) 

with the error formula 

E k = M(h)-P k (h) 

= II M (MlI( A -*•■)*(*)• ( 21 ) 

i=0 i=0 

The next lemma gives an interpolation series for the resolvent using the original k + 1 
interpolation points and setting all of the higher-order interpolation points equal to zero. 
For convenience we shall use the notation M(0) = Mo. Note that if all of the interpolation 
points are set equal to zero the analysis would be that of the Taylor series. 

Lemma 2 Let ho* . . M hk be given and set h m = 0 for all m > k. For 

IM < “in Jj'w - A|, (22) 

a£A{A) 

M may be expanded as 


+oo m m-1 

M(h) = Y. HT n«MlI(‘- M- 

m= 0 i=0 i=0 

Proof. Let l > k. By (21), 

M{h)~ ^(-jrn^ol jv-hi) 


m= 0 


t=0 


i=0 


= (-rf»n«(wn(* -*<)»(*) 

i=0 i= 0 

= (-i) Hi n*wn(‘-w‘) n *•<* 

1=0 i=0 »=fc+ 1 

= { (-»' +1 n ^ (M Tk h - M)M(h)} (hM 0 y- k . 

I i=0 i=0 ) 


(23) 
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But ( hM 0 ) e ~ k -*• 0 as l -* +oo if and only if p(hM 0 ) < 1, which is the well-known 
convergence requirement for a geometric series. From the definition of Mo, we have Mo = 
( jul - A ) -1 . Hence, 


p(jhMo) 


\h\/ min 
xeA(A) 

\h\/r . 


\ju - \\ 


where 


r = min I ju — A|. 
A€A(A) 


(24) 

□ 


This lemma is also important in the development of a pole and zero cancelling routine. 

4 Pole and Zero Cancellation 

Polynomial interpolation approximation works well unless the LTI system being analyzed 
has poles or zeros near the imaginary axis. Such poles and zeros are called resonant poles 
and resonant zeros. The following examples provide the general idea of the effect. 

Example: The deleterious effect of poles and zeros can be illustrated by means of a scalar 
rational function example. Consider 

/(*) = 1 + 2 * + 3 *!i = 1 + 3x + 6x 2 + . . . (25) 

1 - x 

We can use a polynomial approximation to evaluate this function at various values of x. 
Suppose that we choose a second-order polynomial approximation: 

/(x) = 1 + 3x -|- 6x 2 . 

If we evaluate / for x = 0.01 and x = 0.99, we get the approximations 

/(0.01) = 1.0306, 
and 

/(0.99) = 9.8506, 

respectively. If we compare these to the actual values, 

/( 0.01) = 1.0306061 
and 

/(0.99) = 592.03, 
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we can see that as we approach a pole, a much higher-order approximation is required in 
order to get even modest accuracy. 

However, if initially we eliminate the pole before we make our calculations for values near 
x = 1, the accuracy of the method increases dramatically. Again, use a second-order 
approximation with pole cancellation, and we get 

/(*) = (1 - = (1 + 2x + 3a: 2 ). 

i. “■ X 

After evaluating /, we let 


As can be seen in this case, the second-order interpolation is exact. In most cases, however, 
only a marked increase in accuracy is realized. 


In order to cancel a pole in our frequency response, we write 


M(h) = 


(, jh + ju - A) 


and then find a polynomial approximation of (jh+jv— X)M(h). Therefore, our interpolation 

_ Go + G\(h — ho) + • • • + Gk( h — ftp) • • •(/* — hk- 1) 

+ tl ) = * 


(26) 


(jh + ju- A) 

where the coefficient matrices are for a system devoid of the resonance problem. The 
following lemma shows how to compute the new coefficient matrices while preserving the 
form of the interpolating series. 


Lemma 3 Let ho , . . . , ft* be given and set h m = 0 for all m > k. For |/»| < r, where r is 
defined in (24), define the coefficient matrices Fm * implicitly via 


n -foo m—1 

m + ju - A t)M{h) = £ F$P n<*- *»■) • 




m=0 i=0 


(27) 


Then 


ifP = (-;)” - < 28 ) 

t=0 

and 


4”' = (i*» + - K) 4"' ,) + i- ffc* 1 . 

where we define J = 0 for all l. 


m = 0, 1, . 


(29) 
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(30) 


Proof. Equation (28) is immediate from (23). By (27), 

+00 771 — 1 +O0 771 — 1 

uh + jw - a.) £ /!"-■> n (*-*.)■ e iK*- w- 

771=0 1=0 771=0 1=0 

The assumptions that h m = 0 for all m > k and |h| < r ensure that the series in (30) are 

absolutely convergent. We may thus rearrange the left-hand summation as follows: 

+OQ 771 — 1 

(jh + ju-X n )^F^ n (h-hi) 

771=0 1=0 

= E W*- + i" ~ x » + i(* - *-)) Jt"' 0 II (* - *f) 

771=0 1=0 

+OQ 771—1 

= E (jh m + ju - A n )4 n_1) n - w 

771=0 1=0 

+ OQ 771 

771=0 1=0 

+oo m— 1 

= E + JU - A n )il n_1) n ( h - M 

771=0 1=0 

+OC 771—1 

+Ej f t' 1 n (*-*<) 

771=0 1=0 

+00 7 71 — 1 

= E ((>*« + - A n)4 n_1) + jF^) n - m- 


771=0 


i=0 


Comparison with (30) gives (29). □ 

If we need to cancel resonant zeros, we then need to find a polynomial approximation of 
(jh+ju-z) • The following lemma illustrates how this is done. 

Lemma 4 Let ho, hi, . . ., h* be given and set h m = 0 for all m> k. For |h| < r, where r 
is defined in ( 24 ), define the coefficient matrices implicitly via 


+0O 


m— 1 


M(h) = J](j'h + ju - z t ) E D m ] II ” /l *)' 


r=i 


771=0 


1=0 


Then 


and 




t=0 


bL*' = (Bi"- 1 ' - iD™, )/(jh + ju - z.), 

where we define ] — 0 for all t. 


m = 0, 1, . . . 


(31) 

(32) 

(33) 
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Proof-. The proof is similar to that of the preceding lemma except that we start with 
the identity 


(jh + ju - An) £ Dtf n> - hi) = £ ^ _1) II(*- *•■) 

m=0 

and continue from there. 


m-l 

n< 

»=0 


+oo 

£ 

m=0 


m— 1 


i=0 


(34) 

□ 


5 Frequency Response Interpolation Algorithm 

Step 1 Solve for Xq in 

(j(h 0 + u)I- A)X 0 = B , (35) 

and then solve recursively for X\ , . . . , Xk in 

(j(h m + u)I - A)X m = - jX m _! . (36) 


Step 2 Let X$ = X m and define 

X$ = (jhm + ju - A,)^" 1 ) + jX %. Zp , 0 < m < k , (37) 

with Xt\ = 0 for 0 < l < n. 

Step 3 Let Xm ^ = Xm^ and define 

X£M = - jX£M)/(jh m + ju-*t) , 0 < m < k , (38) 

with X l _\ = 0 for 0 < t < l. 

Step 4 Form the coefficient matrices Go, . . G* via 

G = CX£M . (39) 


Step 5 

G(u + h) = (Go + Gi(h - ho) H h G*(/i — ho) • • -(h - i))) 

nLi(j“+j*-»<) no 

mu. 

Remark 

The method used to solve the recursive linear systems in the first step of the algorithm 
depends on the initial structure of the LTI system being investigated. If the system has an 
exploitable structure such as sparsity, an algorithm that exploits that particular structure 
will be used. If no such structure exists, an initial similarity transformation, most likely to 
upper Hessenberg form, will be applied to the system. 
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6 Interpolation Point Selection 

The placement of the interpolation points is of great importance in getting a good ap- 
proximation to the frequency response. We have tested three simple methods to place the 
interpolation points: linear, loglinear, and Chebyshev. We have also tested placement using 
the a priori information of the pole locations. 

Since frequency response is usually plotted against frequency on a log scale, the use of lin- 
early spaced interpolation points does not usually perform well. It places too many points 
at the end of an interval. Both the loglinear and the Chebyshev interpolation point place- 
ments have shown promise. The loglinear placement technique usually gives an excellent 
approximation in the beginning to the middle of an interval, but sometimes can fail mis- 
erably at the end of an interval. The Chebyshev interpolation points (see [8]) spread the 
approximation error fairly evenly across the interval. However, several times the error of 
the Chebychev selection, although acceptable, is larger than that of the acceptable range 
of a loglinear interpolation of the same size. Currently, we are investigating possible hybrid 
techniques to exploit the best of both placement schemes. 

In the cases where we have tried placing interpolation points with the knowledge of the poles 
and zeros of the system the results have been mixed in comparison to the two previously 
mentioned techniques. What has been learned is that under no circumstances should the 
interpolation points be the same as the resonant frequency of a resonant pole or zero. How- 
ever, placing an interpolation point near the resonant frequency improves the approximation 
significantly. 

7 Conclusion 

In this paper we have presented a rational interpolation method for computing the frequency 
response of a system. A significant computational savings can be achieved over several of 
the current methods for computing a frequency response. An error analysis for the method, 
together with other details, can be found in [3]. 

The method presented in this paper avoids the numerical problem of subtraction of near 
equal quantities in the difference terms by using the resolvent identity of Theorem 1. Also, 
simple pole and zero cancellation techniques significantly increase the accuracy of the algo- 
rithm. 

We are currently writing a software package to implement the algorithm in this paper. In 
addition, we are extending this algorithm for use with descriptor systems. 
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